A First-Principles Approach to Insulators in Finite Electric Fields 
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We describe a metliod for computing the response of an insulator to a static, homogeneous electric 
field. It consists of iteratively minimizing an electric enthalpy functional expressed in terms of 
occupied Bloch-like states on a uniform grid of k points. The functional has equivalent local minima 
below a critical field £c that depends inversely on the density of k points; the disappearance of the 
minima at £c signals the onset of Zener breakdown. We illustrate the procedure by computing the 
piezoelectric and nonlinear dielectric susceptibility tensors of III-V semiconductors. 



0^ 



C/3 



I 

o 
o 



> 
cn 
cn 

in 
o 

(N 
O 



I 

o 
o 



X 
S3 



PACS numbers: PACS: 77.22. Ch, 78.20.Bh, 42.70.Mp 

The response of insulators and semiconductors to ex- 
ternal electric fields is of fundamental as well as practical 
interest. It determines their dielectric, piezoelectric, and 
ferroelectric behavior. Much current technological inter- 
est is focused on the use of static fields to tune prop- 
erties such as the dielectric function in the RF and mi- 
crowave region or the index of refraction in the optical 
region. Although a sophisticated physical understand- 
ing of electric field effects has emerged ^, ^, for the 
most part this has not translated into tractable compu- 
tational schemes applicable to periodic solids. Efficient 
first-principles methods for computing derivatives of the 
total energy of solids with respect to a macroscopic field £ 
at £ = do exist ||, §|, 0] . While for many applications 
such perturbation approaches are adequate, their exten- 
sion to nonlinear order is awkward, and for some studies 
(e.g., field- induced structural phase transitions (|]) it is 
essential to perform calculations directly at finite fields. 

The main difficulty is the nature of the scalar potential 
"— 5-r", which is nonperiodic and unbounded from below. 
The fact that it destroys the periodicity of the crystal po- 
tential means that methods based on Bloch's theorem do 
not apply. Moreover, as a result of the unbounded na- 
ture of the perturbation, the energy can always be low- 
ered by transferring charge from the valence states in one 
region to conduction states in a distant region. This is 
the intrinsic dielectric breakdown caused by interband 
(or Zener) tunneling |[ ^ ^ . In many practical situ- 
ations Zener tunneling is negligible on the relevant time 
scale, and for relatively small fields the system remains in 
a polarized long-lived resonant state. This is the state we 
would like to describe. However, the absence of a well- 
defined ground state invalidates the variational princi- 
ple that underlies the usual time-independent electronic- 
structure methods and leads to problematic "runaway" 
solutions in implementations of such approaches pO| . 

The few first-principles methods that have been pro- 
posed for dealing with finite fields in solids have had lim- 
ited success. The supercell "sawtooth" approach be- 
comes prohibitively expensive for all but the simplest sys- 
tems. A significant advance, rooted on the "modern the- 
ory of polarization" ||l^ , was the development of a real- 
space method based on truncated field-polarized Wan- 
nier functions [[l3[ , which removed the need for supercells; 



however, the full first-principles implementation [ p4[ was 
hindered by convergence problems and proved too cum- 
bersome to find widespread use. 

In this Letter we propose an alternative variational ap- 
proach. It is based on the minimization of an electric en- 
thalpy functional J- with respect to Bloch-like functions, 
where J- is comprised of the usual Kohn-Sham energy 
ii^KS and a field coupling term — Pmac ■ £■ Here Eks and 
the polarization Pmac are expressed in terms of a set of 
field-polarized Bloch functions, the latter via the Berry- 
phase theory of polarization . Although for £ ^ the 
Bloch functions arc not eigenstates, they form an appro- 
priate representation of the electronic system. We justify 
this approach, showing that a suitable choice of Brillouin 
zone sampling prevents runaway solutions. We demon- 
strate its implementation into a standard electronic band 
structure program by computing piezoelectric and lin- 
ear and nonlinear dielectric properties of III-V semicon- 
ductors, finding good agreement with experiment. Our 
method opens up new possibilities for first-principles in- 
vestigation of electric-field effects in condensed matter. 

We solve for field-polarized Bloch functions tpnki'r) = 
e*'' '' u„k(r) [where Unk(r) — Mnk(r -I- R)] by minimizing 
the electric enthalpy functional introduced in Ref. m , 



T[u„i^; £] = £^KS - Pi] 



(1) 



where is the primitive cell volume and Pmac = Pion + 
Pel is the macroscopic polarization. In a continuous-A: 
formulation, Pd is --/e/(27r)'^ times the sum of valence- 
band Berry phases iQ J^^dk{unk\i'^ii\unk) (/ is the 
spin degeneracy and e > 0). However, as we show below, 
it is essential to use a formulation in terms of a mesh 
of Nk — Ni X N2 X k points in the Brillouin zone 
(BZ). Then Eks = (f/Nk) Enj("nk, | i?Ks(kj)|M„k,) and 
Pel ■ bi = ife/n) ip'^l' where 
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Imln Yl det S(k^\]<fl^ (2) 
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is the string-averaged discretized Berry phase along the 
direction of primitive reciprocal lattice vector [ p^ . 
Here 5„m(k, k') = (u„k|umk'), n and m run over the 
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M occupied bands, and N^'f^ ~ N2 x N3 is the number 
of strings along bi, each containing A^i points k^"'^' = 

k|^^ + (j7^i)bi. 

It is implicit in this formulation that when £ ^ 
the electronic structure can continue to be represented 
in terms of field-polarized Bloch-like functions V'nki even 
though they are no longer eigenstates of the Hamiltonian. 
It is well known that when describing an occupied sub- 
space, one has the freedom to carry out an arbitrary 
unitary transformation among the states used to repre- 
sent it. In this spirit, we assume that the density matrix 
can be written as p{r,r') = (1/iVfe) X;„k Ck(i"')V'nk(r), 
where n runs over the same number M of Bloch-like 
states at all k. Then p{r,r') = p(r -I- R, r' -I- R) and 
it follows that the charge density and other local quan- 
tities are periodic under translation by a lattice vector 
R. These are familiar properties of an insulating ground 
state, and they turn out to hold for the field-polarized 
state as well, since: (i) If one starts with an insulat- 
ing ground state and applies a homogeneous electric field 
with arbitrary time dependence, the occupied manifold 
preserves the above "insulating-like" properties at later 
times, (ii) A state that minimizes ^ is a stationary solu- 
tion to the time-dependent Schrodinger equation in the 
presence of a static field; since it can be obtained by adi- 
abatically turning on the field, it is guaranteed by (i) 
to have the above "insulating-like" properties. Proofs of 
(i-ii) are not difficult; details will be given elsewhere 

Usually, the transition to the discrete k space is intro- 
duced for merely computational reasons. Here, on the 
contrary, the discrete k formulation is required to elimi- 
nate the possibility of runaway solutions, i.e., to allow for 
stable stationary solutions to exist, unaffected by Zener 
charge-leakage. This is consistent with previous work 
showing that in the continuous-fc limit there are no sta- 
tionary solutions to the static electric field problem (for 
an overview, see Sec. II. D of Ref. [^). To understand 
why the discretization procedure endows J- with minima, 
it is useful to think of it as "bending" space into a finite 
ring: a uniform mesh of Ni x N2 x N3 k points amounts 
to imposing periodic boundary conditions - which have a 
ring topology - over a supercell of dimensions Li = NiOi 
(i=l,2,3). For a given /c-point mesh, J- will have min- 
ima only if £ is small enough that Zener tunneling is 
effectively suppressed. This should happen as long as 
the distance across which the electrons would have to 
tunnel in order to lower their energy is larger than the 
ring perimeter Li = Niai. By thinking of the field as 
spatially tilting the energy bands, one arrives at the con- 
dition £ ■ Sii < £c • Sbi where e£c • ~ i?gap/-^i and Sgap 
is a representative energy gap. We have confirmed this 
behavior in a one-dimensional three-band tight-binding 
model jl^ by studying the stability of the field-polarized 
solutions and by checking that, for a given k point mesh, 
Egiip/eaNk is a reasonably good estimate of £c- 

The polarized state below £c has additional insulating- 
like properties, namely the absence of a steady-state cur- 



rent and the localization of the Wannier functions to a 
small portion of the ring |]l6[ |l^ . This state is related 
to the zero-field ground state by a continuous "deforma- 
tion." Such "polarized manifolds" have been discussed 
in the literature [Q, ||, || for infinite crystals, where they 
were characterized as long-lived resonances. Instead, for 
our finite ring the state obtained by minimizing J- is truly 
stationary, as discussed above. 

By virtue of the nature of the Berry phase term, the 
functional J- cannot be recast as the expectation value of 
a Hermitian operator. Because that term contains over- 
laps between the states at neighboring k points, even 
in a tight-binding model without charge self-consistency, 
the problem must be solved self-consistently among all k 
points. This breakdown of Bloch's theorem is the price 
to pay for handling, within periodic boundary conditions, 
a field whose scalar potential breaks translational invari- 
ance. Indeed, the Berry-phase term in is the proper re- 
placement of the usual scalar potential term e£ ■ (r) when 
using a ring topology. (Alternatively, one can switch to a 
vector potential formalism, which restores translational 
invariance to the Hamiltonian at the expense of rendering 
the static field problem time-dependent Q.) 

Let us now describe the minimization algorithm. We 
have chosen an iterative "band-by-band" conjugate- 
gradients method [ p^ in which each occupied state M„k 
is updated in sequence, although other schemes may be 
equally suitable. The many-electron state of interest vi- 
olates inversion symmetry but not time-reversal sym- 
metry; the latter can be used, together with any £- 
preserving point-group operations, in reducing the BZ. 
The few differences with respect to a normal ground-state 
calculation stem from the — fiPmac • £ term, as follows. 
The gradient |G„k) = S!F/S{unk\ becomes 

|G„k) = -f ^Ks(k)|u„k) + E ^ X (3) 

E [k„,k(;>)^™\(k,k^^) - |u„^,^«)5-i(k,k«)_ , 

■m— 1 

where k^-* = k ± /Ni and use was made of Eq. (88) of 
Ref. . By standard manipulations this is converted 
into a preconditioned conjugate-gradients search direc- 
tion |Dnk) orthonormalized to the occupied manifold at 
k. The trial updated state is \u„]s}{0) = cos^?|7i„k} + 
sm6\Dnk). We search J-{0) for a minimum, replace |ii„k} 
by |S„k), and go on to the next band. 

However, the behavior of ^{6) is unconventional, as is 
illustrated in Fig. ||. The -Eks(^) contribution (dashed 
line in Fig. ^ is the usual one; it is periodic with period 
TT and has an amplitude proportional to E^l/Nk- How- 
ever, —V,P^ac{d)-£ (dotted line) has a secular component 

arising from the fact that Pmac changes by l/N^f^ times 
a "quantum of polarization" AP = feR/Q where R 
is a (usually nonzero) lattice vector, as 9 —>■ O + n. To un- 
derstand this, consider one phase (3 = Imlndet ^(k, k') 
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strain, it follows from the expression 




FIG. 1: The electric enthalpy T (solid line) and its compo- 
nents -Eks (dashed line) and — fiPmac-f (dotted line), plotted 
as a function of the parameter Q that controls the update of 
a polarized Bloch state in a conjugate-gradients step. 



contributing to Eq. (|^). Its 6'-dependence is /3(0) — 
Im In {A cos 6* -I- B sin 6') , where A = det S" and B = AqXS 
are complex constants and S is obtained from S by re- 
placing the n-th row with (I?,ik|wmk')- It is then easy to 
see that /3(0) progresses by ±7r as Q increases by tt (the 
sign depending on the sign of ImA*_B). In total there 
are six such contributions for each k; for £ along b^, 
the two involving neighboring k' along the string direc- 
tion i contribute to an average slope —jeS ■ a.i/Nf'T: of 
— J7Pinac(^) • 5 as a function of 9, as shown by the dotted 
line in Fig. ^ If f is not too large, we simply carry out the 
update by stepping to the nearest local minimum of T{9) . 
(Local minima separated by tt are equivalent.) However, 
J- {9) loses its minima when € gets too big. This occurs 
for e£ ■ a,; ~ E^g^p/Ni, which is precisely the heuristic 
condition for the onset of Zener tunneling on a ring. 

One possible concern with the present method is that 
it imposes a minimum mesh spacing that can be used 
for a given £. This difficulty can be circumvented in 
practice by decomposing the fine k mesh into a set of 
sparser submeshes, computing the Berry-phase terms on 
each submesh, and then averaging over all of them. Of 
course, -Eks can be computed on the fine mesh. 

We now turn to the computation of forces and stresses 
at £ ^ 0. According to the HcUmann-Feynman argu- 
ment at a stationary point of T the force Fj = 
—dT/dYj becomes simply —dT/drj, i.e., the implicit de- 
pendence via the wave functions can be dropped. Eq. (|^) 
has no such explicit dependence, and so does not con- 
tribute to Fj. Thus, aside from the trivial ionic core con- 
tribution eZj the force is given by the standard £—Q 
HcUmann-Feynman expression arising from Eks alone. 

As for the macroscopic stress, similar arguments yield 
(^afj = {^l^)dJ- / drjafi^ where ry is a homogeneous strain. 
However, the electric boundary conditions under which 
the strain derivative is taken must be specified carefully. 
When the cell is deformed according to ^ (1 -I- 'q)sLi, 
we can hold fixed either the macroscopic field £, or the 
potential drop across each lattice vector, Vi = £ ■ a.;. 
As the Berry phases ip^^^ do not depend explicitly on the 



27rPei • £ = ^ (£ • a,)(Pei • b^) = ^ 5] (£ • a,) ^ 



i=l 



(4) 



that 9(r2Poi • £)/driaf3 = when V is fixed (the same 
holds true for the ionic term, which can also be recast in 



terms of a phase ip. 



I). As a result, the stresses in 



the two cases are related by 



'al3 



^(V) 



^^f„(a.),(^lf (5) 
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SO that the pressures differ by (Pmac • ^)/3. The stress 
cr(^) is given in terms of the polarized Bloch states by the 
same expression as the stress in a zero-field ground-state 
calculation; it is a symmetric (torque-free) bulk quantity. 
On the contrary, cr^^^ generally has an antisymmetric 
part, and moreover it depends on the choice of branch 
cut when evaluating the multivalued bulk polarization. 
(In the context of a finite crystallite, the torque and ex- 
tra stress in cr^^^ can be regarded as arising from forces 
exerted on the polarization- induced surface charges by a 
field that is held fixed in the "laboratory frame" |£l|.) It 

is straightforward to show that 0)^^^ = da'gj /d£a is the 
so-called "improper" piezoelectric tensor pO|, whereas 
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dV, 



(6) 



is the "proper" tensor. 

The scheme outlined above was first validated on a one- 
dimensional tight-binding model where the results 
were found to agree with the results of the Wannier-based 
approach It was then implemented in ABINIT [ p2[ , 
a fully self-consistent pseudopotential code. To illustrate 
the utility of the method, we have calculated by finite- 
differences dielectric and piezoelectric constants of sev- 
eral III-V semiconductors. That is, we increase £ in 
small increments and compute the changes in the re- 
sulting forces, stresses, and polarizations, with internal 
displacements and strain either clamped or undamped 
as appropriate. The Born effective charge is eZ*^^ = 
dFjp / d£a- (Contrary to previous finite-difference ap- 
proaches, we compute it as the derivative of a force with 
respect to £, not polarization with respect to displace- 
ment.) The dielectric constant is = ^ap + Xap^ where 
Xa0 = (1/eo) d{Pmi,c)a/d£i3 and cq is the vacuum permit- 
tivity. If the ions are kept fixed, this yields the electronic 
contribution 60©; if both electrons and ions are allowed 
to relax in response to the field, the static dielectric con- 
stant estat is obtained. The quadratic susceptibility is 
Xal-y = i2/eQ)d^{P^i,c)a/d£f3d£j, and we have computed 
it keeping the ions fixed. In the zinc blende structure, the 
only nonzero independent components of these tensors 
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GaAs 


AlAs 


GaP 


MP 


a (a.u.) 


10.45 


10.59 


10.11 


10.24 




(10.68) 


(10.69) 


(10.28) 


(10.33) 


■'-'cation 


2.00 


2.14 


2.10 


2.24 




(2.07) 


(2.18) 


(2.04) 


(2.28) 


Coo 


11.9 


9.6 


9.4 


8.1 




(10.9) 


(8.2) 


(9.0) 


(7.5) 


^static 


13.5 


11.5 


11.2 


10.2 




(13.2) 


(10.1) 


(11.1) 


(9.8) 


X<^) (pm/V) 


134 


64 


66 


39 




(166) 




(74) 




(aVe)7i4 


-0.40 


-0.10 


-0.25 


0.05 




(-0.32) 




(-0.18) 




-1.42 


-1.40 


-1.35 


-1.31 



TABLE I: Computed dielectric and piezoelectric properties of 
III-V semiconductors. Parentheses denote experimental data 
quoted in Refs. @, |[ ||, ||]. 

are Zli, en, and Xi23- The "proper" piezoelectric coeffi- 
cient c{^^ is computed using Eq. (^), with both clamped 

ilii) ^ii'i undamped (714) ions. 

The calculations were performed at the theoretical lat- 
tice constant a using an energy cutoff of 10 Ha. We 
checked that our k point sampling was very well con- 
verged at 16 X 16 X 16 points in the full BZ. With this mesh 
spacing we find critical fields of the order of lO^V/cm, 
and the finite-difference field step size was approximately 
1/10 of this value. We checked that our values for Z*, 



7i4, and 7^4 essentially coincide, for any given mesh of 
k points, with those computed using the approach of 
Refs. |2^. We have also computed Z* and eoo by 
treating the electric field via density-functional perturba- 
tion theory (DFPT) Q. In the limit of a dense mesh the 
two approaches yield the same results; the discrepancies 
that occur for sparser meshes can be attributed to the 
different ways in which derivatives with respect to k are 
handled in each case P| . Our results for the piezoelectric 
coefficients and for x a-re also in good agreement with 
experiment and with previous calculations using different 
methods |, |, 

All of the quantities reported in Table I could have 
been obtained using DFPT methods. However, a consid- 
erable gain in convenience is afforded by computing them 
using simple finite differences of £. For example, the cal- 
culation of x*^^^ by DFPT is quite tedious and requires 
a special-purpose program, while x*^"-* of arbitrary order 
are easily computed using the present approach. 

To summarize, we have presented a practical first- 
principles scheme for computing the electronic structure 
of insulators under a finite dc bias. The algorithm is ide- 
ally suited for implementation in a standard electronic 
structure code and its computational cost is comparable. 
Dielectric polarization, forces and stresses are straightfor- 
wardly obtained as byproducts of the calculation. The 
extension of this approach to time-dependent fields will 
be discussed in a future communication. 

This work was supported by NSF Grant DMR- 
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